clear all
clc
close all

load AALabData
Row=row;
Column=column;

load AANewData
Row=[Row;row'];
Column=[Column;column'];

load AAMTData
Row=[Row;row];
Column=[Column;column];

%% Make data

ROW=zeros(486,9);
for i=1:486
    ROW(i,:)=reshape(Row{i},[9,1])';
end
[~,I]=unique(ROW,'rows');

Row=Row(I,:);
Column=Column(I,:);
n=length(Row);

%% Parameter grid

lambda=0:0.001:0.01;
u=length(lambda);

%% Generate synthetic data

Q=1000;
fakedistr=cell(1,Q);
ErrorUniform=zeros(1,Q);
ErrorLogit=zeros(1,Q);
complete=9999*ones(1,Q);
fitlambda=zeros(1,Q);

load FM

for q=1:Q
    % Generate fake distributions
    fakedistr{q}=genfake(Row,Dominant,Dominated);
    Freq=fakedistr{q};
        
        % Predict uniform
            
        err=zeros(1,n);
        for i=1:n
            err(i)=Freq(i,1)*log((Freq(i,1)+eps)/(1/3))+Freq(i,2)*log((Freq(i,2)+eps)/(1/3))+Freq(i,3)*log((Freq(i,3)+eps)/(1/3));
        end
        ErrorUniform(q)=mean(err);
        
        % Run PCHM     
        
        levR{1}=1/3*ones(1,3);
        levC{1}=1/3*ones(1,3);
        
        lambdaerr=zeros(1,u);
        for k=1:u
            lambdaerr(k)=logitlev1fixedlambda(lambda(k),Row,Freq);
        end
        fitlambda(q)=find(lambdaerr==min(lambdaerr));
        ErrorLogit(q)=min(lambdaerr);
        
end

bestlambda=lambda(fitlambda);

'Restrictiveness'
r = mean(ErrorLogit)/mean(ErrorUniform)

'Standard Error'
sampleVarLL1 = mean((ErrorLogit-mean(ErrorLogit)).^2);
sampleVarUniform = mean((ErrorUniform-mean(ErrorUniform)).^2);
covar = mean((ErrorLogit-mean(ErrorLogit)).*(ErrorUniform-mean(ErrorUniform)));

SE = sqrt((sampleVarLL1 - 2 * r * covar + r^2 * sampleVarUniform) / mean(ErrorUniform)^2) * (1/sqrt(Q))
